Energy Diffusion and Heat Conduction Near Thermal Equlibrium 
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We derive analytically a dynamical equality relating the mean square displacement of energy diffusion near 
equilibrium to the autocorrelation function of total heat flux fluctuation in equilibrium. This relation is a di- 
rect consequence of the microscopic energy conservation law. We numerically verify the dynamical equality 
for typical nonlinear lattice systems, and also show some special consequences of this dynamical equality, in 
particular, the connection between anomalous heat conduction and anomalous energy diffusion without a priori 
assumption for the underlying diffusion processes. 
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It is well known that normal heat conduction is related to 
normal diffusion. For example, in a ID continuous system, 
the Fourier's law of heat conduction states that j{x,t) = 
—KdxT{x,t), where j{x,t) is the local heat flux density, 
T{x,t) the local temperature and k the thermal conductivity. 
From energy continuity equation dte{x,t) + dxj{x,t) — 0, 
one obtains dte{x,t) — KdxT{x,t). Assuming local ther- 
mal equilibrium exists and temperature change is small, the 
local energy density can be simply related to the local temper- 
ature as e{x,t) = cvT{x,t), where cy is the specific heat. 
Therefore, one reaches the familiar heat equation describing 
normal diffusion dtT{x,t) = DdxT{x,t), where the diffu- 
sion constant D is linked with heat conductivity k through 
D = k/cv[1]. 

The discovery of anomalous heat conduction [2] in ID 
Fermi-Pasta-Ulam (FPU) and FPU-aUke lattices has attracted 
enormous theoretical [3-6] and experimental [7, 8] studies 
of the length-dependent heat conductivity in low dimensional 
systems. In such systems, the asymptotic heat conductivity 
can be expressed as 



(1) 



On the other hand, for the diffusion processes of particles 
[9], which could be heat carriers, it is established that the 
mean square displacement (MSD) of particles increases with 
time as 



x^{t) oc f^. 



(2) 



Since the Fourier's law (a — 0) can be connected to nor- 
mal diffusion (7 = 1), it is commonly believed that anoma- 
lous heat conduction (a ^ 0) should somehow be related to 
anomalous diffusion (7 ^ 1). The pioneering works in this 
direction have given rise to insightful albeit controversial an- 
alytical results [10, II] for billiard models. Further numerical 
investigations of energy diffusion in ID lattice systems [12- 
15] confirmed the relation a = 7 — 1 between heat conduction 
and diffusion as predicted by Ref . [11]. However, the analyt- 
ical derivation of this relation is only available for the simple 



non-interacting billiard model under a strong assumption that 
the moving particles obey a pre-defined Levy walk distribu- 
tion. Therefore, this result is still restricted to dealing with 
the particle diffusion cases. Actually, similar relation between 
electric conduction and diffusion has also been investigated in 
the context of the diffusion of charged particles [16]. 

However, energy transport in discrete lattices is a com- 
pletely different phenomenon. The heat carriers are phonons 
which are defined in momentum space and do not move in the 
real space [17]. Therefore the definition of "displacement" of 
energy is not clear. We do not have the quantity "mean square 
displacement" of energy using which we are able to character- 
ize energy diffusion in these lattice cases. As a consequence, 
the previous inspiring works that tried to bridge energy diffu- 
sion and heat conduction from the aspect of particle diffusion 
have only limited appUcation. A fundamental theory for lat- 
tice models in this area is still lacking, although some heuristic 
results are available [13, 15]. 

Therefore, to build a general and rigorous theory connect- 
ing energy diffusion and heat conduction that does not rely on 
any a priori assumption on the details of the system and can 
be applied to any model as well as lattices is still challenging 
and remains to be one of the most outstanding problems in 
statistical physics. 

In this Letter, based on a rigorous linear response descrip- 
tion of energy diffusion, we derive analytically a dynamical 
equality connecting the MSD of the nonequilibrium energy 
diffusion x^ (t) and the autocorrelation function of total heat 
flux, C(t), in equiUbrium: 
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2C{t) 
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where cy is the heat capacity per unit length, ks the Boltz- 
mann constant, and T the temperature. Our derivation does 
not rely on any special details of the model or any assump- 
tion about the diffusion process. It is a direct consequence of 
microscopic energy conservation. In this sense, the equality 
we derived is a fundamental and intrinsic relation between 
energy diffusion and heat conduction. We further numeri- 
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cally verify our dynamical equality by using two paradigmatic 
ID nonlinear lattice systems, i.e., 0^ lattice showing normal 
heat conduction and quartic FPU-/? (q'FPU-/?) lattice exhibit- 
ing anomalous heat conduction. 

Derivation For simplicity, we limit ourselves to a ID con- 
tinuous system. The generalization of the derivation to higher 
dimensional and discrete lattice cases is straightforward. In 
thermal equilibrium, we have the microscopic energy conti- 
nuity equation: 



dte{x,t)+d^j{x,t)^0. 



(4) 



where j{x,t) denotes the local heat flux density at position 
X at time t. We have used the energy density fluctuation 
e{x, t) = Ae(a;, t) = e{x, t) — {e{x, t)) to replace the energy 
density e(x,t) in the above equation because the difference 
between them is a time-independent constant. (•) denotes the 
canonical ensemble average. 

The continuity equation (4) connects local energy density 
fluctuation with local heat flux density. To find the further 
connection between correlation functions of these two quanti- 
ties, we multiply Eq. (4) by e{x' , t') and j{x\ t') respectively, 
and take the ensemble average: 

dt {e{x, t)€{x', t')) + d^{jix, t)eix', t'}) = 0, (5) 
dt{e{x,t)j{x',t'))+d^{j(x,t)j{x',t'))^0. (6) 

For a homogeneous system in thermal equilibrium, the 
correlation functions {e{x,t)e{x' ,t')), {j{x,t)e{x' ,t')), 
{e{x, t)j{x' ,t')) and {j{x, t)j{x' , t')) only depend on spatial 
separation and time differences. Therefore we can define 
Cab{x-x' ,t-t') =CBA{x'-x,t'-t) = {A{x,t)B{x\t')), 
where A, B takes the form of e, j. It is then easy to verify 
that the cross correlation terms Cje{x,t) and Cej{x,t) are 
invariant under the transformation of {x,t) {~x,—t). 
Therefore, we have Cje {x — x' ,t — t') = C^j {x — x' ,t — t'). 
By performing dt to Eq. (5) and dx to Eq. (6), and subtracting 
the cross terms, we obtain: 



dfcUx,t) = dlCjj{x,t). 



(7) 



As we shall see soon, the rh.s. of Eq. (7) is actually con- 
nected to the equilibrium autocorrelation of total heat flux 
C{t) and the l.h.s. is closely connected to the MSD x^{t) 
of the nonequilibrium energy diffusion. 

The autocorrelation function of total heat flux is defined as 



Cit)= lim i(J^(f)Ji(0)) 



(8) 
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where Jl = /_^/2 ^(•'^i t^Ax is the total heat flux for ID sys- 



tem with length L. Thus Cit) can be expressed as 

r.L/2 .L/2 
-L/2 J -L/2 



1 rL/2 .L/2 

C(t) = lim - / / Cjj{x - x',t)dxdx' 



(9) 



where J is short for throughout the paper Now, the 
connection between the rh.s of Eq. (7) and C{t) has been 
established. 

Before we come to the connection between the l.h.s of Eq. 
(7) and the MSD x^{t), we need to obtain the energy distri- 
bution p{x, t) describing nonequilibrium energy diffusion so 
that the MSD, or the variance of the energy distribution, can 
be defined as 



x^{t) 



x'^p{x, t)dx. 



(10) 



It should be stressed that the MSD x-^ (t) is not measured in 
the canonical ensemble as denoted by (.). Instead, it is defined 
in a non-stationary irreversible diffusion process described by 
the probability distribution of energy p{x, t). 

The scenario of energy diffusion is given as follows: firstly, 
an energy perturbation is imposed on the initially equilibrated 
system at temperature T. Then after removing the perturba- 
tion, the energy diffusion process is just the relaxation of the 
energy excitation. Following Ref. [4], the external perturba- 
tive Hamiltonian is introduced as 



Hextit) = -ei-t) / dxe{x,t) 



^ AT{x) 
T '' 



(11) 



under the assumption of local thermal equilibrium, where 6{t) 
is the Heaviside function which indicates the perturbation is 
removed at time t = after it has been applied from the in- 
finite past. As we will see soon, AT (a;) describes the energy 
profile before relaxation. The form of AT{x) will not affect 
our main result, namely, Eq. (3) is generic and does not de- 
pend on the specific initial energy profile. 

Provided |AT(a:;)| is sufficiently small, we can apply the 
Hnear response theory [18] to study the energy distribution 
during the diffusion process such that for t > the response 
of e{x,t) to the perturbation is 



5e{x,t)^[ e{~t')dt' I dx'^ee[x~x' ,t^t') 

-^^"di'l dx'<i>,,ix-x',t')ATix'). 



(12) 



^es{x,t) is known as the response function which relates 
to the equilibrium correlation of energy fluctuation as [18] 

(^^,{x,t) = -/3t (e(0,0)e(x,t)) = -PTCee{x,t) , where ■ 
denotes time derivative and (5t = {ksT)^^ . Using the condi- 
tion C^^{x, oo) — 0, we obtain 



(Se(x,t)= j dx'Cee{x 



AT{x') 



(13) 



which describes the evolution of energy profile in the diffusion 
process. 

At time i = 0, the total energy excited in the system 6E = 
J Se{x, 0)dx reads 



SE =JJ dxdx'CUx~x',0)^^ = J 



cvAT{x')dx' , 
(14) 
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where / Cee{x, 0)dx = ksT'^cy is used [19]. This result im- 
plies an expected result that the total energy excited is the 
accumulation of local energy perturbation. Actually, when 
the interaction in the system is short ranged, the simulta- 
neous correlation of energy density can be approximated as 
Cee(a;, 0) ^ d{x), then the energy profile at t = can be 
simplified as de{x,0) = cv^T{x). In this aspect, Ar(a;) 
describes the initial excitation profile of energy diffusion. 

In view of the analysis above, the normalized energy distri- 
bution in the diffusion process is written as 



p{x,t) 



Se{x, t) 
5E 



^ I Ax'C, 



{x- x',t)ATix'), (15) 



where TV = kBT^cy J AT{x')dx' is the normalization con- 
stant. From energy conservation, it is expected that A{t) — 
J p{x, t)dx should always be unity. Using Eq. (7), this can 
be proved as follows. Since d'^A{t)/dt'^ ~ JJ d'^Cjj{x — 
x',t)AT{x')dxdx' /J\f EE 0, together with dA{0)/dt = 
from the even parity of Cee{x, t) with respect to t and A{0) = 
1 by the definition of p{x,t), we get dA{t)/dt = and 
Ait) = 1. 

For a special perturbation AT(x) oc S{x), Eq. (15) 
reduces to p{x,t) = C^e{x,t)/kBT'^cv- This result is 
reminiscent of a heuristic definition for discrete lattices 
p{i,t) = {AHo{0)AH,{t)) / {AHo{0)AHom [13] pro- 
vided (Ai7o(0)Ai7o(0)) = ksT'^cv, where H^{t) is the lo- 
cal energy of site i in the lattice, the discrete analogue of 
e{x,t). However, this condition implies (A_ffo(0)AiJi(0)} = 
for i 7^ or its continuous version C^e{x,0) = 
kBT^cvS{x), which has more restrictions on the system, e.g., 
short ranged interaction and special definition for the local en- 
ergy [21]. Therefore, our result Eq. (15) is more general and 
rigorous. 

To describe the evolution of the MSD, from Eq. (7), (10) 
and (15), we obtain 



dV(t) 

di2 



■i // x'&iC,,(x 



x' ,t)AT{x')dxdx' . (16) 



Because J Cjj [x, t)dx exists as indicated in Eq. (9), we 
have the boundary conditions limj._).oo x^dxCjj {x, t) = and 
limaj^oo xCjj{x, t) = 0. By performing integration by parts 
twice on Eq. (16) with respect to x, we finally arrive at our 
central reiM/f Eq. (3) for arbitrary AT{x). 

This dynamical equality is essentially the equation of mo- 
tion for the MSD of energy diffusion. Its initial conditions are 
specified as ^(0) = x^Cee{x - x' ,Q)AT{x')dxdx' /N 
and dx'^ (0) /dt = 0. The first one reduces to x'^ (0) = if 
Cee(a;, 0) and AT{x) are both ^-functions. The second one 
results from the even parity of C^eix, t) with respect to t pro- 
vided it exists [22]. These initial conditions imply that the 
lowest two terms of Taylor series of x"^ (t) at i = should 
be ao + a2t^, with oq = ^^(O) and 02 = C{0)/kBT^cv. 
Therefore, realistic energy diffusion process should start as a 
ballistic transport. Similar phenomenon for particle diffusion 
has already been observed in experiments recently [23]. 



As the central result of our work, Eq. (3) relates the MSD 

(t) of the nonequilibrium energy diffusion and the equilib- 
rium autocorrelation function of total heat flux, C{t). This 
relation bridges the heat conduction and energy diffusion in 
the sense that: (1) heat conduction is related to the autocorre- 
lation function C{t) [24]; (2) energy diffusion processes can 
be classified through the distinct behaviors of MSD x^{t) [9]. 

We should stress that Eq. (3) is derived without any explicit 
prior assumption on the details of the system and diffusion 
process. The derivation is independent of the initial condition 
of energy diffusion. It is an intrinsic relation purely coming 
from the microscopic energy conservation law. 

Simulation In order to verify Eq (3), we perform equilib- 
rium numerical simulations on ID nonlinear lattices with a 
general dimensionless Hamiltonian: 



H 



Pi 
2 



+ Viq^-q,-,) + Uiq^), (17) 



where and Pi denotes the displacement and momen- 
tum, V{q) the inter-particle potential, U{q) the on-site 
potential and Hi the local Hamiltonian or energy. We 
adopt the definition in Ref. [4] for the local heat flux 
ji{t) ~ —qidV{qi — qi-i) /dqi and the total heat flux 
J{t) = The energy distribution is calculated ac- 

cording to the discrete form of Eq. (15), i.e., p{i,t) = 
{AHoiO)AH,{t)) /kBT\vyithAH,it) = H,{t)^{H,{t))- 
The MSD is calculated as x'^{t) = Ei*^P(«:0- The Boltz- 
mann constant fcs is set as unit. We note that different def- 
initions of local energy and heat flux do not affect the veri- 
fication as long as they are defined consistently according to 
energy conservation law. The Hamiltonian is integrated us- 
ing the fourth order symplectic scheme CSABA2 [25, 26] with 
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FIG. 1: Numerical verification of Eq. (3) for (a) the diffusive 
lattice (N = 501, T = 0.441, cv = 0.885 [21]); and (b) the su- 
perdiffusive q¥?\]-l3 lattice (N = 601, T = 0.020, cv = 0.750). 
ft = 0.1 is used for small t while ft = 5 (gFPU-/3) or 10 (0") for 
large t in order to reduce the fluctuation. 
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time step At = 0.05 which usually conserves the total en- 
ergy with a relative deviation of the order of 10^® after a total 
time 10^. The second derivative of the MSD is calculated as 
<Pl^{t)/de = \^[t + h)- 2^{t)+li^[t - h)]/h^ where 
h is the time interval. 

Our numerical verification is performed on the 0'' lattice 
with Viq) = U{q) = 9^/4 and the qFPV-/3 lattice 

with V{q) — (7^/4, U{q) = 0, respectively. The former is 
a paradigmatic lattice system showing normal heat conduc- 
tion [27, 28]. The latter is just the high temperature limit 
of the FPU-/3 model, a well studied lattice system exhibit- 
ing anomalous heat conduction [29-32]. The q'FPU-/? lattice 
has temperature independent heat capacity cy = 3/4 which 
can be easily deduced from the generalized equipartition theo- 
rem (zdzH) = ksT for any canonical coordinate z = qi,Pi', 
while the heat capacity for (f)^ lattice is temperature depen- 
dent which is calculated from cv{T) = dT{E{T)) /N , where 
^^(r) is the total energy of the system at temperature T [21]. 
The numerical results in Fig. (1) show perfect agreement with 
the dynamical equality (3) for both models, which demon- 
strate its universality. 

Application In the following, we will show some special 
corollaries of Eq. (3), with only the help of Green-Kubo (GK) 
formula [24]. 

(1): If the energy diffusion is normal, we have 
limt_j.oo a;^(t)/2t = D, where D is a finite constant called 
thermal diffusivity. From the GK formula and noticing the 
initial conditions for x'^{t), we obtain the thermal conductiv- 
ity: 



.normal 



C{t) cv dj;2(t) 



We shall emphasize that although Eq. (18) is well accepted, 
it was obtained by the aids of a presumed empirical Fourier's 
law, or equivalently the heat equation, as shown at the begin- 
ning of this paper However, we have demonstrated here that 
it is a rigorous result from Eq. (3) and does not rely on the 
Fourier's law any more. 

Particularly, if we assume an exponential relaxation of the 
flux autocorrelation function C{t) = Cqc^*/'^ similar to the 
(/)^ lattice in Fig. (la), with r the characteristic relaxation 
time, the MSD can be obtained by integrating Eq. (3) twice, 
reading 



[t] = 2D t - r ( 1 - c 



(19) 



The asymptotic behavior is x'^(t) ^ 2Dt with diffusivity 
D = Cqt/ (kBT'^cv). The thermal conductivity is k = 
C(t)dt/{kBT^) = CoT/ikBT'^), so that k and D satisfy 
the corollary Eq. (18). We also notice that in the short time 
limit t <^ T, Eq. (19) implies ballistic transport x'^{t) — v^t^, 
where v = ^ Co/iksT'^cv) can be regarded as the speed of 
the collective energy propagation in the ballistic limit. The 
heat conductivity can be expressed as k = cyu^r which is 
reminiscent of the kinetic theory. It can be noticed that x'^ (t) 



in Eq. (19) has exactly the same time dependence as the MSD 
of Brownian motion described by the Langevin equation [18]. 

(2): If the energy diffusion is superdijfusive, then asymptot- 
ically we have x'^{t) ~ f with 1 < 7 < 2. From the equality 
Eq. (3) we get C{t) - with 6 = - 2. The GK formula 
is not applicable due to the divergence of Jp°° C{t)dt. Thus 
we get anomalous heat conduction with infinite conductivity 
in the thermodynamic limit. However, if one is interested in 
a length-dependent heat conductivity, one can follow the heat 
conduction review literatures [4, 5] to put a cut-off time to the 
GK formula. Although not rigorously proved, it is commonly 
reasoned as that for a finite system, sound waves travels to the 
boundaries at a constant speed c,, which will lead to a fast 
decay of correlations at time t ^ L/cg [5]. In this case, our 
central result Eq. (3) gives 



^supcr 



kuT^ 



, (20) 



t~L/c 



which leads to the asymptotic behavior of the length- 
dependent thermal conductivity k^"*""^ ^ L" with 



1 = 7-1. 



(21) 



This reproduces the relation a = 7 — 1 derived for a billiard 
model with the particles follow the a priori Levy walk process 
[11]. On the contrary, our derivation is a direct consequence of 
the microscopic energy conservation law and does not rely on 
any assumption on the specific models and the diffusion pro- 
cesses. It is fulfilled in the whole time range of the diffusion 
process. In this sense, our dynamical equation (3) is intrinsic 
and more fundamental than a simple relation between 7 and 
a. 

(3): Let us look into a situation in which the MSD x'^{t), 
based on the equilibrium energy spatiotemperal correlation 
Eq. (15), has the asymptotic form x^{t) ^ V with < 7 < 1, 
which corresponds to the subdijfusive case. From Eq. (3), we 
get C{t) - 7(7 - VjV-'^. The exponent 5 = 7 - 2 < -1 so 
that C [t) is integrable over [0, 00) and the GK formula is still 
applicable. From Eq. (3), we have 



..sub 
K 



lim 



cy dx^jt) 

2 dt 



lim t'' 



0. 



(22) 



This is an expected result which indicates a subdiffusive sys- 
tem is a perfect thermal insulator in the thermodynamic limit. 
If we again adopt the "cut-off-time" reasoning for finite sys- 
tem, Eq. (20) and (21) will still be valid and the heat conduc- 
tivity will decay as, k™^ L''^^, asymptotically. 

It is interesting to note that when t is large, the autocorre- 
lation function C{t) should be negative since its coefficient 
^ 7(7 — 1) < 0. This is in stark contrast to the case where 
C{t) ~ with (5 < — 1 but has a positive coefficient, which 
will result in a finite conductivity and thus implies normal dif- 
fusion [4, 33]. 

In conclusion, we have analytically derived an exact dy- 
namical equality, Eq. (3), for all isotropic systems which 
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bridges energy diffusion and heat conduction in terms of the 
MSD of the nonequiUbrium energy diffusion, x'^{t), and the 
equihbrium autocorrelation function of total heat flux, C{t). 

This relation is a natural consequence of microscopic en- 
ergy conservation law and is thus intrinsic and more funda- 
mental than the previously derived relation (1 = 7 — !. In 
principle, it can be applied to isotropic systems in any dimen- 
sion 

From this equality, we are able to recover the theories for 
normal and anomalous heat conduction. The validity of this 
intrinsic relation has been verified with two paradigmatic 1 D 
nonlinear lattice systems, i.e., the 0^ lattice which shows nor- 
mal heat conduction and the qFPV-fi lattice which exhibits 
anomalous heat conduction. The perfect agreements between 
theory and numerical simulations demonstrate its universality. 

At last, we would like to mention that the present results are 
derived in the linear response framework. Whether there ex- 
ist universal relations beyond linear response regime to unify 
heat conduction and energy diffusion is still an open question 
and requires further investigation. 
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